      MODULE WATER_DENSITY
      use PHYS_CONSTANTS, only : &
     & row0
      implicit none

!     The coefficients of UNESCO formula
!     for water density dependence on temperature

      real(8), parameter:: a0 = 800.969d-7
      real(8), parameter:: a1 = 588.194d-7
      real(8), parameter:: a2 = 811.465d-8
      real(8), parameter:: a3 = 476.600d-10

!     The coefficients of formula for
!     water density dependence on temperature and salinity
!     (McCutcheon et al.,Water Quality and Maidment. 
!      Handbook on hydrology, 1993)
      
      real(8), parameter:: A11 = 288.9414d0
      real(8), parameter:: A12 = 508929.2d0
      real(8), parameter:: A13 = 68.12963d0
      real(8), parameter:: A14 = 3.9863d0
      real(8), parameter:: alpha1 = 8.2449d-1
      real(8), parameter:: alpha2 = 4.0899d-3
      real(8), parameter:: alpha3 = 7.6438d-5
      real(8), parameter:: alpha4 = 8.2467d-7
      real(8), parameter:: alpha5 = 5.3675d-9
      real(8), parameter:: beta1 = 5.724d-3
      real(8), parameter:: beta2 = 1.0227d-4
      real(8), parameter:: beta3 = 1.6546d-6
      real(8), parameter:: gamma1 = 4.8314d-4
        
      contains
      REAL(8) FUNCTION WATER_DENS_T_UNESCO(Temp)

!     The function WATER_DENS_T_UNESCO
!     returns the density of water, kg/m**3, 
!     as a function of temperature 
!     according to UNESCO formula
      
      implicit none

      real(8), intent(in):: Temp
      
      WATER_DENS_T_UNESCO = &
     &  row0*(1+a0+a1*Temp-a2*Temp**2+a3*Temp**3)
      END FUNCTION WATER_DENS_T_UNESCO


      REAL(8) FUNCTION WATER_DENS_TS(Temp,Sal)

!     The function WATER_DENS_TS
!     resturns the density of water, kg/m**3,
!     as a function of temperature and salinity
!     according to
!     (McCutcheon et al.,Water Quality and Maidment. 
!      Handbook on Hydrology, 1993)

!     Input variables:
 
!     Temp --- water temperature, deg C
!     Sal  --- water salinity,    kg/kg
      
      implicit none

      real(8), intent(in):: Temp
      real(8), intent(in):: Sal

      real(8) Sal_g_kg,A,B,C

!     Converting salinity units from kg/kg to g/kg
      Sal_g_kg = Sal*1.d+3

!     The first term, dependent on temperature      
      WATER_DENS_TS = row0 * &
     & ( 1.-(Temp+A11)*(Temp-A14)**2/(A12*(Temp+A13) ) )

      A =  alpha1         - alpha2*Temp    + alpha3*Temp**2 &
     &    -alpha4*Temp**3 + alpha5*Temp**4

      B = -beta1          + beta2*Temp     - beta3*Temp**2

      C = gamma1

!     Adding the second term, dependent on temperature and salinity
      WATER_DENS_TS = WATER_DENS_TS + &
     & A*Sal_g_kg + B*Sal_g_kg**1.5 + C*Sal_g_kg**2.

      END FUNCTION WATER_DENS_TS


      REAL(8) FUNCTION DDENS_DTEMP(Temp,Sal)

!     The function DDENS_DTEMP
!     returns the derivative of water density
!     on temperature, kg/(m**3*C)
      
!     Input variables:
 
!     Temp --- water temperature, deg C
!     Sal  --- water salinity,    kg/kg

      implicit none

      real(8), intent(in):: Temp
      real(8), intent(in):: Sal      

!     Note: the terms with salinity are omitted,
!           i.e. this is the derivavive for fresh water

      DDENS_DTEMP = &
     & -(row0/A12)*(Temp-A14)/( (Temp+A13)**2 ) * &
     & ( (Temp+A13)*( (Temp-A14)+2.*(Temp+A11) ) - &
     &  (Temp+A11)*(Temp-A14) )

      END FUNCTION DDENS_DTEMP
      

      REAL(8) FUNCTION DDENS_DSAL(Temp,Sal)

!     The function DDENS_DSAL
!     returns the derivative of water density
!     on salinity, , kg/(m**3*kg/kg)
      
!     Input variables:
 
!     Temp --- water temperature, deg C
!     Sal  --- water salinity,    kg/kg

      implicit none

      real(8), intent(in):: Temp
      real(8), intent(in):: Sal

      real(8) Sal_g_kg,A,B,C

!     Note: the dependency of A and B on temperature
!           is omitted

      A =  alpha1
      B = -beta1
      C =  gamma1
	     
!     Converting salinity units from kg/kg to g/kg
      Sal_g_kg = Sal*1.d+3

      DDENS_DSAL = &
     & A + 1.5*B*Sal_g_kg**0.5 + 2.*C*Sal_g_kg

      DDENS_DSAL = DDENS_DSAL*1.d+3

      END FUNCTION DDENS_DSAL
      
      
      FUNCTION DENS_HOSTETLER(temp)
      
!     Function DENS_HOSTETLER calculates the density of water
!     dependent on temperature only, following Hostetler and Bartlein, 1990

      implicit none

      real(8) :: DENS_HOSTETLER
      
      real(8), intent(in) :: temp
      
!     Parameters      
      real(8), parameter :: row0 = 1.d+3
      real(8), parameter :: temp0 = 3.85
      real(8), parameter :: const = 1.9549d-5
      real(8), parameter :: const_power = 1.68
      
      DENS_HOSTETLER = row0*(1.-const*dabs(temp-temp0)**const_power)
      
      END FUNCTION DENS_HOSTETLER
      END MODULE WATER_DENSITY


      MODULE PHYS_FUNC
      contains
      REAL(8) FUNCTION TURB_DENS_FLUX(tempflux,salflux,Temp,Sal)

!     Function TURB_DENS_FLUX            _____
!     returns the turbulent density flux w'ro' in water
!     at given temperature, 
!              salinity,         ____
!              temperature flux  w'T'
!                                ____
!          and salinity    flux  w's'

!     Input variables:
!     tempflux --- kinematic heat flux, m*C/s
!     salflux  --- salinity       flux, m*(kg/kg)/s
!     Temp     --- temperature        , C
!     Sal      --- salinity           , kg/kg

      use water_density, only : &
     & ddens_dtemp, &
     & ddens_dsal
      
      implicit none
      
      real(8), intent(in):: tempflux
      real(8), intent(in):: salflux
      real(8), intent(in):: Temp
      real(8), intent(in):: Sal

      TURB_DENS_FLUX = &
     & ddens_dtemp(Temp,Sal)*tempflux + &
     & ddens_dsal (Temp,Sal)*salflux

      END FUNCTION TURB_DENS_FLUX

   
!      REAL(8) FUNCTION dirdif()
!      use atmos, only:
!     & shortwave
!     implicit none
!      real(8) cloud
!      real(8) dirdif0,b,S0,sinh0
!      common /cloud/ cloud
!     data cloud /0./

!      SAVE

!      b=1./3.
!      S0=1367.
!      cloud = 0.
    
!      dirdif0 = (shortwave-b*S0*sinh0())/(b*(S0*sinh0()-shortwave))
!      dirdif = dirdif0*(1.-sngl(cloud))
!      dirdif = dmax1(dirdif,0.d0)
      
!      END FUNCTION dirdif


      FUNCTION SINH0(year,month,day,hour,phi)

!     SINH0 is sine of solar angle ( = cosine of zenith angle)   
      
      implicit none
      
      real(8) :: sinh0
      
      integer(4), intent(in) :: year
	  integer(4), intent(in) :: month
	  integer(4), intent(in) :: day
	
	  real(8)   , intent(in) :: hour
	  real(8)   , intent(in) :: phi

      real(8) delta
	  real(8) theta
	  real(8) pi
	  real(8) phi_rad
	
      integer(4) nday

      integer(4), external:: JULIAN_DAY
    
      pi=4.*datan(1.d0)

      nday  = JULIAN_DAY(year,month,day)

      delta = 23.5d0*pi/180.d0*dCOS(2*pi*(float(nday)-173.d0)/365.d0)
      theta = pi*(hour-12.d0)/12.d0
      phi_rad = phi*pi/180.d0
      sinh0 = dSIN(phi_rad)*dSIN(delta) + &
      & dCOS(phi_rad)*dCOS(delta)*dCOS(theta)
      sinh0=dmax1(sinh0,0.d0) 

      END FUNCTION SINH0


      REAL(8) FUNCTION QS(phase,t,p)
      
!     QS - specific humidity, kg/kg, for saturated water vapour

      implicit none
      real(8) t,p,aw,bw,ai,bi,a,b,es
      integer(4) phase
!phase = 1 - liquid, 0 - ice
      aw = 7.6326    
      bw = 241.9    
      ai = 9.5  
      bi = 265.5 
      a  = phase*aw+(1-phase)*ai
      b  = phase*bw+(1-phase)*bi
      es=610.7*10.**(a*t/(b+t))
      QS=0.622*es/p
      END FUNCTION


      REAL(8) FUNCTION ES(phase,t,p)
      
!     ES is pressure of saturated water vapour, Pa    

      implicit none
      real(8) t,p,aw,bw,ai,bi,a,b
      integer(4) phase
!phase = 1 - liquid, 0 - ice
      aw = 7.6326    
      bw = 241.9    
      ai = 9.5  
      bi = 265.5 
      a  = phase*aw+(1-phase)*ai
      b  = phase*bw+(1-phase)*bi
      ES = 610.7*10.**(a*t/(b+t))
      END FUNCTION

      
      REAL(8) FUNCTION MELTPNT(C)
      implicit none
      real(8) C, dtdc
      
      SAVE
      
      dtdc = 66.7
      
      MELTPNT = 0. - C*dtdc
      END FUNCTION MELTPNT
      
      
      FUNCTION WATER_FREEZE_MELT(temperature, grid_size, &
     & melting_point, switch)
     
!     The function WATER_FREEZE_MELT cheks, if the
!     heat storage of a grid cell is larger or equal
!     the latent heat, necessary for phase transition
      
      use PHYS_CONSTANTS, only : &
     & cw_m_row0, &
     & ci_m_roi, &
     & row0_m_Lwi, &
     & roi_m_Lwi
      use NUMERIC_PARAMS, only : &
     & min_ice_thick, &
     & min_water_thick, &
     & T_phase_threshold
     
      implicit none
      
!     Input variables
      real(8), intent(in) :: temperature ! Temperature, Celsius
      real(8), intent(in) :: grid_size
      real(8), intent(in) :: melting_point
      
      integer(4), intent(in) :: switch ! +1 for water -> ice transition
                                       ! -1 for ice -> water transition
      logical :: WATER_FREEZE_MELT
      
      if (switch == +1) then
        if ( ( melting_point - T_phase_threshold - temperature) * &
     &   cw_m_row0*grid_size > min_ice_thick*roi_m_Lwi .or. &
     &   ( melting_point - T_phase_threshold - temperature) > 0.2d0) &
     &   then
          WATER_FREEZE_MELT = .true.
        else
          WATER_FREEZE_MELT = .false.
        endif
      elseif (switch == -1) then
        if ( (temperature - T_phase_threshold - melting_point) * &
     &   ci_m_roi*grid_size > min_water_thick*row0_m_Lwi .or. &
     &   (temperature - T_phase_threshold - melting_point) > 0.2d0) &
     &   then
          WATER_FREEZE_MELT = .true.
        else
          WATER_FREEZE_MELT = .false.
        endif     
      else
        write(*,*) 'Wrong switch in WATER_FREEZE_MELT: STOP'
        STOP
      endif
      
      END FUNCTION WATER_FREEZE_MELT


      SUBROUTINE TURB_SCALES(k_turb_T_flux,T_massflux,eflux0_kinem,h1, &
      & turb_density_flux, Buoyancy0, H_mixed_layer, w_conv_scale, &
      & T_conv_scale)
     
      use ARRAYS, only : &
      & dzeta_05int
      use DRIVING_PARAMS, only: &
      & M
      use PHYS_CONSTANTS, only: &
      & g, &
      & row0

      implicit none

!     Input variables
      real(8), intent(in) :: k_turb_T_flux(:)
      real(8), intent(in) :: T_massflux(:)
      real(8), intent(in) :: eflux0_kinem
      real(8), intent(in) :: h1
      real(8), intent(in) :: turb_density_flux
      real(8), intent(in) :: Buoyancy0
	
!     Output variables
      real(8), intent(out) :: H_mixed_layer
      real(8), intent(out) :: w_conv_scale
      real(8), intent(out) :: T_conv_scale
      
!     Extrnal functions
      real(8), external:: DZETA
      
!     Help variable
      integer i, maxlocat      
	
      maxlocat = 1
      do i = 2, M
	    if (k_turb_T_flux(i) + T_massflux(i) > &
        & k_turb_T_flux(maxlocat) + T_massflux(maxlocat)) maxlocat = i
      enddo
      
      H_mixed_layer = dzeta_05int(maxlocat)*h1

      w_conv_scale = (Buoyancy0*H_mixed_layer)**(1.d0/3.d0) 

      T_conv_scale = -eflux0_kinem/w_conv_scale
     
      END SUBROUTINE TURB_SCALES
	

      REAL(8) FUNCTION W_SEDIM()
!     The function W_SEDIM calculates the sedimentation speed
!     of hydrosol particles      
      implicit none
	
!     Currently the sedimentation speed is assumed to be constant,
!     however it is dependent on many properties of both the
!     particles and the liquid media	
      W_SEDIM = 1.0d-5
	
      END FUNCTION W_SEDIM


      FUNCTION WATER_ALBEDO(sinh0)
      implicit none

      real(8) :: WATER_ALBEDO

!     Input variables
      real(8), intent(in) :: sinh0

!     Local variables
      real(8), parameter :: const1 = 0.05d0
      real(8), parameter :: const2 = 0.15d0

      WATER_ALBEDO = const1/(sinh0 + const2)

      END FUNCTION WATER_ALBEDO

      
	  FUNCTION EXTINCT_SNOW(snow_density)
	  implicit none

      real(8) :: EXTINCT_SNOW

!     Input variables
      real(8), intent(in) :: snow_density ! Snow density, kg/m**3

!     Local variables
      real(8), parameter :: const1 = 0.13d0
      real(8), parameter :: const2 = 3.4d0
      real(8), parameter :: extinct_snow_max = 1.d+7

      EXTINCT_SNOW = dexp(-(const1*snow_density+const2) )
      
!     For Sparkling2002-2005 experiment
!     EXTINCT_SNOW = dexp(-extinct_snow_max)


 	  END FUNCTION EXTINCT_SNOW


      REAL(8) FUNCTION UNFRWAT(T,i)
      use driving_params
      use arrays
      
!    CALCULATION OF LIQUID WATER CONTENT IN FREEZING SOIL
!    T - DEG. C; WLM0,WLM7,WLMAX - KG/KG

      implicit none
      real(8) T
      integer(4) i

      unfrwat = (WLM0(i)-WLM7(i))*dexp(T/3.) + WLM7(i)
      !unfrwat = 0 
      
      RETURN
      END FUNCTION UNFRWAT


      REAL(8) FUNCTION WL_MAX(por,rosdry,wi)
      use PHYS_CONSTANTS, only: &
      & row0, &
      & row0_d_roi
      implicit none
      real(8), intent(in):: por,rosdry,wi

      WL_MAX = por*row0/(rosdry*(1 - por)) - wi*row0_d_roi
      END FUNCTION WL_MAX


      REAL(8) FUNCTION WI_MAX(por,rosdry)
      use PHYS_CONSTANTS, only: &
     & roi
      implicit none
      real(8), intent(in):: por,rosdry

      WI_MAX=por*roi/(rosdry*(1-por))
      END FUNCTION WI_MAX

      
      FUNCTION SOIL_COND_JOHANSEN(wl,wi,rosdry,por)
      
      use PHYS_CONSTANTS, only: &
      & row0, roi, &
      & row0_d_roi, roi_d_row0, &
      & lamw0, lami
      
      implicit none
      
      real(8) :: SOIL_COND_JOHANSEN
      
!     Input variables
      real(8), intent(in) :: wl, wi ! liquid water and ice content 
                                    ! (in respect to dry soil mass), kg/kg
      real(8), intent(in) :: rosdry ! dry soil (soil particles) density, kg/m**3
      real(8), intent(in) :: por ! soil porosity, m**3/m**3
      
!     Local variables
      real(8) :: water_vol_ratio, ice_vol_ratio ! water and ice volume ratio 
                                                ! (in respect to bulk soil volume), m**3/m**3
      real(8) :: waterice_vol_ratio ! total volume ratio of water and ice, m**3/m**3
      real(8) :: Kersten ! Kersten number, n/d
      real(8) :: CK_const, CK_const1, CK_const2
      real(8) :: lambda_sat, lambda_dry ! heat conduction coefficients for saturated
                                        ! and dry soil, W/(m*K)
                                        
      real(8) :: water_sat_ratio, ice_sat_ratio                       
                       
      real(8), save :: CK_consts(1:4) ! Cote and Konrad (2005) constants, n/d
      real(8), save :: CK_consts1(1:3) ! Cote and Konrad (2005) constants, W/(m*K)
      real(8), save :: CK_consts2(1:3) ! Cote and Konrad (2005) constants, n/d
      real(8), save :: quartz_ratio ! quartz content of the total solids content
      real(8), save :: lambda_quartz ! heat conduction coefficient for quartz, W/(m*K)
      real(8), save :: lambda_othmin ! heat conduction coefficient for non-quartz
                                     ! minerals, W/(m*K)
      real(8), save :: lambda_solids ! heat conduction coefficient for soild part
                                     ! of the soil, W/(m*K)
      
      logical, save :: firstcall = .true.
      
      if (firstcall) then
        CK_consts(1) = 4.6d0 ! for gravel and coarse sand
        CK_consts(2) = 3.55d0 ! for medium and fine sand
        CK_consts(3) = 1.9d0 ! silty and clay soils
        CK_consts(4) = 0.6d0 ! organic fibrous soils
        
        CK_consts1(1) = 1.7d0 ! for crashed rocks
        CK_consts1(2) = 0.75d0 ! for mineral soils
        CK_consts1(3) = 0.3d0 ! organic fibrous soils
        
        CK_consts2(1) = 1.8d0 ! for crashed rocks
        CK_consts2(2) = 1.2d0 ! for mineral soils
        CK_consts2(3) = 0.87d0 ! organic fibrous soils
        
        quartz_ratio = 0.1d0
        lambda_quartz = 7.7d0
        if (quartz_ratio > 0.2d0) then
          lambda_othmin = 2.d0
        else
          lambda_othmin = 3.d0
        endif
        lambda_solids = lambda_quartz**quartz_ratio * &
        &               lambda_othmin**(1.d0-quartz_ratio)
      endif
      
!     Convertation from mass ratios to volume ratios      
      water_vol_ratio = wl / &
      & (por*(wl + row0/roi*wi + row0/rosdry))
      ice_vol_ratio = wi / &
      & (por*(wi + roi/row0*wl + roi/rosdry))
      
      waterice_vol_ratio = water_vol_ratio + ice_vol_ratio
      
      CK_const = CK_consts(3) ! silty and clay soils are assumed
      Kersten = CK_const*waterice_vol_ratio / &
      & (1.d0 + (CK_const - 1.d0)*waterice_vol_ratio)
      
      water_sat_ratio = por*water_vol_ratio/waterice_vol_ratio
      ice_sat_ratio = por*ice_vol_ratio/waterice_vol_ratio
      
!     This is the original formula from Johansen (1975) extended for the case
!     with the ice content
      lambda_sat = lambda_solids**(1.d0 - por) * &
      & lamw0**(water_sat_ratio)*lami**(ice_sat_ratio)
      
      CK_const1 = CK_consts1(2) ! mineral soils are assumed
      CK_const2 = CK_consts2(2) ! mineral soils are assumed
!     Cote and Konrad (2005) formula for heat conduction coefficient of dry soil
      lambda_dry = CK_const1*10.d0**(-CK_const2*por)
      
!     The normalized soil conductivity concept by Johansen (1975)      
      SOIL_COND_JOHANSEN = (lambda_sat - lambda_dry)*Kersten + lambda_dry
      
      if (firstcall) firstcall = .false.
      END FUNCTION SOIL_COND_JOHANSEN


      FUNCTION REACPOT_ARRHEN &
      & (delta_E, temp, temp0, reacpot_arrhen0)

!     The FUNCTION REACPOT_ARRHEN calculates the reaction potential
!     according to Arrhenius equation

      use PHYS_CONSTANTS, only : &
      & R_univ

      implicit none      
      
!     Input variables
      real(8), intent(in) :: delta_E ! activation energy, J/mol
      real(8), intent(in) :: temp    ! temperature, Kelvin
      real(8), intent(in) :: temp0   ! reference temperature, Kelvin
      real(8), intent(in) :: reacpot_arrhen0 ! reaction potential at the
                                             ! reference temeperature, mol/(m**3*s)
                                             
!     Output variables
      real(8) :: REACPOT_ARRHEN
      
!     Local variables and constants
      
      REACPOT_ARRHEN = &
      & reacpot_arrhen0*dexp(delta_E/R_univ*(1./temp0 - 1./temp))
      
      END FUNCTION REACPOT_ARRHEN
      
      
      FUNCTION HENRY_CONST(henry_const0, temp_dep, temp_ref, temp)
      
      ! Function HENRY_CONST calculates the Henry constant of a substance for a given temperature
      
      implicit none
      
      real(8) :: HENRY_CONST
      
      ! Input variables
      real(8), intent(in) :: henry_const0 ! Henry constant at the reference temperature, mol(m**3*Pa)
      real(8), intent(in) :: temp_dep ! Temperature dependence (enthalpy solution devided by universal gas constant), K
      real(8), intent(in) :: temp_ref ! Reference temperature, K
      real(8), intent(in) :: temp ! Temperature, K
      
      HENRY_CONST = henry_const0*dexp(-temp_dep*(1./temp-1./temp_ref))
      
      END FUNCTION HENRY_CONST
      
      
      FUNCTION DIFF_WATER_METHANE(temp_C)
      
      ! Function DIFF_WATER_METHANE calculates the molecular diffusivity 
      ! of methane dissolved in liquid water, m**2/s
      ! (Broecker and Peng, 1974)
      
      implicit none
      
      real(8) :: DIFF_WATER_METHANE
      
      ! Input variables
      real(8), intent(in) :: temp_C ! temperature, degrees Celsius

      ! Local variables
      real(8), parameter :: const1 = 9.798d-10
      real(8), parameter :: const2 = 2.986d-11
      real(8), parameter :: const3 = 4.381d-13
      
      DIFF_WATER_METHANE = const1 + const2*temp_C + const3*temp_C*temp_C
       
      END FUNCTION DIFF_WATER_METHANE
      
      
      FUNCTION DIFF_AIR_METHANE(temp_C)
      
      ! Function DIFF_WATER_METHANE calculates the molecular diffusivity 
      ! of methane dissolved in liquid water, m**2/s
      ! (Lerman, 1979)
      
      implicit none
      
      real(8) :: DIFF_AIR_METHANE
      
      ! Input variables
      real(8), intent(in) :: temp_C ! temperature, degrees Celsius

      ! Local variables
      real(8), parameter :: const1 = 1.875d-5
      real(8), parameter :: const2 = 1.3d-7
      
      DIFF_AIR_METHANE = const1 + const2*temp_C
       
      END FUNCTION DIFF_AIR_METHANE

      
      FUNCTION DIFF_WATER_CARBDI(temp_C)
      
      ! Function DIFF_WATER_CARBDI calculates the molecular diffusivity 
      ! of crabon dioxide dissolved in liquid water, m**2/s
      ! (Broecker and Peng, 1974)
      
      implicit none
      
      real(8) :: DIFF_WATER_CARBDI
      
      ! Input variables
      real(8), intent(in) :: temp_C ! temperature, degrees Celsius

      ! Local variables
      real(8), parameter :: const1 = 9.39d-10
      real(8), parameter :: const2 = 2.671d-11
      real(8), parameter :: const3 = 4.095d-13
      
      DIFF_WATER_CARBDI = const1 + const2*temp_C + const3*temp_C*temp_C
       
      END FUNCTION DIFF_WATER_CARBDI

      
      FUNCTION DIFF_AIR_CARBDI(temp_C)
      
      ! Function DIFF_WATER_CARBDI calculates the molecular diffusivity 
      ! of carbon dioxide dissolved in liquid water, m**2/s
      ! (Lerman, 1979)
      
      implicit none
      
      real(8) :: DIFF_AIR_CARBDI
      
      ! Input variables
      real(8), intent(in) :: temp_C ! temperature, degrees Celsius

      ! Local variables
      real(8), parameter :: const1 = 1.35d-5
      real(8), parameter :: const2 = 0.9d-7
            
      DIFF_AIR_CARBDI = const1 + const2*temp_C
       
      END FUNCTION DIFF_AIR_CARBDI
      
      
      FUNCTION SCHMIDT_NUMBER_METHANE(tempC)
      
      implicit none
      
      real(8) :: SCHMIDT_NUMBER_METHANE
      
      ! Input variables      
      real(8), intent(in) :: tempC
      
      ! Local variables
      real(8), parameter :: const1 = 1.898d+3
      real(8), parameter :: const2 = -1.101d+2
      real(8), parameter :: const3 = 2.834
      real(8), parameter :: const4 = -2.791d-2
      
      SCHMIDT_NUMBER_METHANE = &
      & const1 + const2*tempC + const3*tempC**2 + const4*tempC**3
      
      END FUNCTION SCHMIDT_NUMBER_METHANE
      
      
      FUNCTION SCHMIDT_NUMBER_OXYGEN(tempC)
      
      implicit none
      
      real(8) :: SCHMIDT_NUMBER_OXYGEN
      
      ! Input variables      
      real(8), intent(in) :: tempC
      
      ! Local variables
      real(8), parameter :: const1 = 1.8006d+3
      real(8), parameter :: const2 = -1.201d+2
      real(8), parameter :: const3 = 3.7818
      real(8), parameter :: const4 = -4.7608d-2
      
      SCHMIDT_NUMBER_OXYGEN = &
      & const1 + const2*tempC + const3*tempC**2 + const4*tempC**3
      
      END FUNCTION SCHMIDT_NUMBER_OXYGEN
      
      
      FUNCTION SCHMIDT_NUMBER_CARBDI(tempC)
      
      implicit none
      
      real(8) :: SCHMIDT_NUMBER_CARBDI
      
      ! Input variables      
      real(8), intent(in) :: tempC
      
      ! Local variables
      real(8), parameter :: const1 = 1.911d+3
      real(8), parameter :: const2 = -1.137d+2
      real(8), parameter :: const3 = 2.967
      real(8), parameter :: const4 = -2.943d-2
      
      SCHMIDT_NUMBER_CARBDI = &
      & const1 + const2*tempC + const3*tempC**2 + const4*tempC**3
      
      END FUNCTION SCHMIDT_NUMBER_CARBDI
      
      
      FUNCTION GAS_WATATM_FLUX &
      & (tempC,wind10,surf_conc,partial_pressure,henry_const,gasindic)
            
      ! Subroutine GAS_WATATM_FLUX calculates the upward gas flux at the water-air interface, mol/(m**2*s)
      ! following formulations by (McGillis et al., 2000; Cole and Caraco, 1998; Riera et al., 1999) 
      ! and others (described in Wania, 2007)
            
      implicit none
      
      real(8) :: GAS_WATATM_FLUX
      
      ! Input variables
      real(8), intent(in) :: tempC ! water surface temperature, Celsius
      real(8), intent(in) :: wind10 ! wind speed at 10 m above the surface, m/s
      real(8), intent(in) :: surf_conc ! gas concentration at the water surface, mol/m**3
      real(8), intent(in) :: partial_pressure ! partial pressure of a gas in the atmosphere, Pa
      real(8), intent(in) :: henry_const ! Henry constant of a gas, mol/(m**3*Pa)
      
      integer(4), intent(in) :: gasindic ! gas indicator: methane - 8, oxygen - 9, carbon dioxide - 10
      
      ! Local constants
      real(8), parameter :: constvel1 = 5.75d-6 ! SI units
      real(8), parameter :: constvel2 = 5.97d-7 ! SI units
      real(8), parameter :: Schmidt_scale = 600.
      real(8), parameter :: wind10_power = 1.7
      
      integer(4), parameter :: methane_indic = 8
      integer(4), parameter :: oxygen_indic = 9
      integer(4), parameter :: carbdi_indic = 10
      
      ! Local variables
      real(8) :: Schmidt_number
      real(8) :: piston_velocity, piston_velocity600
      real(8) :: surf_conc_atmequil
      
      piston_velocity600 = constvel1 + constvel2*wind10**wind10_power
      if (gasindic == methane_indic) then
        Schmidt_number = SCHMIDT_NUMBER_METHANE(tempC)
      elseif (gasindic == oxygen_indic) then
        Schmidt_number = SCHMIDT_NUMBER_OXYGEN(tempC)
      elseif (gasindic == carbdi_indic) then
        Schmidt_number = SCHMIDT_NUMBER_CARBDI(tempC)
      endif
      piston_velocity = piston_velocity600*dsqrt(Schmidt_scale/Schmidt_number)
      surf_conc_atmequil = partial_pressure*henry_const
      GAS_WATATM_FLUX = piston_velocity*(surf_conc-surf_conc_atmequil)
      
      END FUNCTION GAS_WATATM_FLUX
      
      
      FUNCTION CHARNOCK_Z0(velfrict)

      use PHYS_CONSTANTS, only : &
      & g, niu_atm
      
      implicit none

      real(8) :: CHARNOCK_Z0

      ! Input variables
      real(8), intent(in) :: velfrict

      ! Local variables
      real(8), parameter :: const1 = 0.111
      real(8), parameter :: const2 = 0.0144
      real(8), parameter :: roughness_min = 1.d-4
      real(8), parameter :: roughness_max = 1.1d-1

      CHARNOCK_Z0 = &
      & dmin1(dmax1(const1*niu_atm/velfrict + const2*velfrict**2/g, &
      & roughness_min),roughness_max)

      END FUNCTION CHARNOCK_Z0


      END MODULE PHYS_FUNC
